IMPROVING FERMI ORBIT DETERMINATION AND PREDICTION IN AN 
UNCERTAIN ATMOSPHERIC DRAG ENVIRONMENT 


Matthew A. Vavrina (1) , Clark P. Newman®, Steven E. Slojkowski (3) , 
and J. Russell Carpenter (4) 

<l) a. i. solutions, Inc., 10001 Derekwood Ln., Suite 215, Lanham, MD 20706, 301-286-5105, 

matthew. vavrina@ai-solutions. com, 

(2> a. i. solutions, Inc., 10001 Derekwood Ln., Suite 215, Lanham, MD 20706, 281-244-0788, 

clark.newman@ai-solutions.com, 

<3> a. i. solutions, Inc., 10001 Derekwood Ln., Suite 215, Lanham, MD 20706, 301-286-6637, 

steven.slojkowski@ai-solutions.com 

<4i NASA Goddard Space Flight Center, 8800 Greenbelt Rd, MS 595.0, Greenbelt, MD 20771, 

301-286- 7526, russell. carpenter@nasa.gov 

Abstract: Orbit determination and prediction of the Fermi Gamma-ray Space Telescope 
trajectory is strongly impacted by the unpredictability and variability of atmospheric density and 
the spacecraft’s ballistic coefficient. Operationally, Global Positioning System point solutions 
are processed with an extended Kalman filter for orbit determination, and predictions are 
generated for conjunction assessment with secondary objects. When these predictions are 
compared to Joint Space Operations Center radar-based solutions, the close approach distance 
between the two predictions can greatly differ ahead of the conjunction. This work explores 
strategies for improving prediction accuracy and helps to explain the prediction disparities. 
Namely, a tuning analysis is performed to determine atmospheric drag modeling and filter 
parameters that can improve orbit determination as well as prediction accuracy. A 45% 
improvement in three-day prediction accuracy is realized by tuning the ballistic coefficient and 
atmospheric density stochastic models, measurement frequency, and other modeling and filter 
parameters. 
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1. Introduction 

Atmospheric drag uncertainty of Earth-orbiting spacecraft at low altitudes poses challenges to 
accurate orbit determination (OD) and prediction. Such is the case for the Fermi Gamma-ray 
Space Telescope orbit, or Fermi, which at an operational altitude of roughly 535 km is strongly 
perturbed by atmospheric drag. The components comprising the drag force model are 
challenging to estimate and are even more difficult to predict due to atmospheric density 
variations caused by unpredictable space weather. Drag modeling inaccuracies are further 
exacerbated by an uncertain spacecraft ballistic coefficient (BC). 

For operational OD of the Fermi spacecraft, point solutions from the spacecraft’s Global 
Positioning System (GPS) receiver are processed by the Fermi owner/operator (O/O) using the 
extended Kalman filter (EKF) and backwards smoothing capabilities in Orbit Determination 
ToolKit (ODTK) [1], Predictive ephemerides are generated for conjunction assessment and are 
compared to predictions produced by the Joint Space Operations Center (JSpOC) using the 
Astrodynamic Support Workstation (ASW) software when a close approach with a second object 
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is anticipated. The JSpOC orbit predictions, however, are based on OD using skin radar 
measurements and are expected to be less accurate than predictions based on an orbit state 
derived from GPS data. GPS measurements are, however, subject to elevation-dependent 
ionospheric path delays, which can reduce OD accuracy [2], Significant differences between the 
predicted miss distance of Fermi with a second orbiting object from the O/O and the JSpOC 
solutions were reported in a number of different conjunction reports from the NASA Robotic 
Conjunction Assessment Risk Analysis (CARA) team. Figure 1 illustrates this behavior in four 
separate cases, depicting the evolution of both the 0/0 (green line) and JSpOC (blue line) 
predicted miss distances as functions of day of year leading up to the time of closest approach 
(TCA), based on daily orbit predictions. Additionally, in several instances, the JSpOC prediction 
appears to predict a more consistent miss distance as the TCA nears despite the use of low- 
quality radar measurements for OD. While the predictions eventually converge on a similar miss 
distance at the TCA, differences as high as 10 km at six days from the TCA are observed. The 
disparity between the predictions leads to difficulty in selecting which solution to trust when 
making a decision on whether to execute a maneuver to mitigate a close approach. This 
discrepancy is the primary motivation for this investigation. 




Figure 1: Evolution of miss distance between Fermi and four secondary objects leading up 
to time of closest approach as predicted by 0/0 ODTK and JSpOC ASW. 

This investigation begins with a tuning analysis to help characterize the observed issues, uncover 
any systemic problems with the 0/0 ODTK operational scenario, and to refine and optimize the 
OD-prediction system. Initially, a broad search is conducted across OD-prediction parameters. 
Filter parameter settings and models are widely varied to identify the parameters that have the 
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highest impact on performance and to detect any strong couplings between parameters. The 
principal tuning strategies are found to be those that are related to atmospheric drag modeling, 
such as the stochastic process model and associated parameters employed for the estimation of 
the BC and atmospheric density. A broad search produces potential leads in terms of improved 
prediction accuracy, and those basins of attraction are investigated in closer detail with a more 
refined and strategic variation of the parameters. The tuning process includes iterating with the 
following strategies: 

1) thinning GPS measurements that are processed in the EKF, 

2) adjusting measurement noise in the filter, 

3) varying the initial drag coefficient of the filter, 

4) evaluating different ballistic coefficient stochastic process models and parameters 
(Gauss-Markov, random walk, and Vasicek models), 

5) modifying ballistic coefficient uncertainty, 

6) inflating atmospheric density correction uncertainty, 

7) investigating different atmospheric models, 

8) testing interpolated versus step-function produced geomagnetic indices. 

The propagation accuracy of the Fermi trajectory is found to be highly sensitive to atmospheric 
drag modeling and uncertainty due to Fermi’s low altitude. Fermi experiences relatively large 
drag perturbations and any mis-modeling of the drag strongly affects prediction accuracy. The 
focus of the tuning is, in turn, determining BC and atmospheric density modeling parameters that 
would best balance definitive OD accuracy with the ability to appropriately maintain the 
estimated value through prediction. Determining the best balance is difficult given the manner in 
which the BC and atmospheric density correction states are modeled as stochastic processes. 
Significant performance improvements are obtained with appropriate tuning. A 45% 
improvement in the mean 3 -day position accuracy and a 26% reduction in the maximum position 
error over the evaluation period are realized by thinning the measurements processed by the 
filter, setting a relatively short half-life for a Gauss-Markov stochastic process for the BC, setting 
a relatively long exponential half-life for the atmospheric density correction stochastic model, 
and inflating the uncertainty of the atmospheric density correction estimated in the filter. 
Additionally, the standard deviation of the maximum 3 -day prediction error decreases by 36%, 
indicating more consistent predictions. 

2. LEO Prediction Accuracy 

Prediction accuracy for any orbit regime is dependent on the models incorporated for estimation 
and propagation. Drag modeling and the associated atmospheric density model are critically 
important for a low-Earth orbiter. Other LEO spacecraft are examined to provide context for 
Fermi orbit prediction accuracy. Drag modeling techniques employed for low-Earth orbiters are 
noted in Sections 2.1 and 2.2, respectively. 

In this study prediction error is measured as the definitive-predictive overlap in position. In the 
absence of truth data, the definitive-predictive overlap provides a reasonable proxy for the error 
from the true orbit state. The magnitude of the difference between the estimated position vector 
and the predicted position vector is referred to as the prediction error. 
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2.1 Hubble Space Telescope Prediction Accuracy 

To appreciate the prediction results from the CARA TCA plots in Figure 1 and to properly 
evaluate the results of any modifications to the operational ODTK scenario, it is important to 
understand the prediction accuracy that should be expected for a LEO satellite at Fermi’s 
altitude. An appropriate case study is the Hubble Space Telescope (HST). While the two 
spacecraft have a different ballistic coefficient, they are in similar circular orbits, where Fermi’s 
perigee altitude is approximately 535 km and HST is slightly higher at a perigee altitude of 
approximately 550 km (note that the altitude for both spacecraft is variable). Accurate orbit 
prediction for HST and Fermi is challenging because of their altitude, where the atmospheric 
drag is difficult to forecast with great accuracy. In general, atmospheric density models are 
typically less accurate the higher the altitude; however, atmospheric uncertainty plays a more 
prominent factor in prediction at lower altitudes where the drag is a major perturbing force. 


The NASA Flight Dynamics Facility (FDF) maintains historical prediction accuracy data for 
HST, serving as a gauge for what might be expected in terms of orbit prediction accuracy for 
Fermi. Historical orbit prediction accuracy for Hubble is also captured in a technical report by 
Nicholson and Dell [3]. The 3-day prediction accuracy for HST from July 1996 through 
December 2012 is plotted in Figure 2 along with the observed FI 0.7 solar flux index. Through 
2011, the FDF performed the OD and prediction using the Goddard Trajectory Determination 
Software (GTDS), which uses batch a least squares approach for orbit estimation. Conspicuous 
in Figure 2 are the relatively high prediction errors during the last solar maximum period 
(indicated by a high F10.7 index) in the early the 2000’s, where there are several 3-day 
prediction accuracies over 20 km and many days over 5 km. 
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Figure 2: 3-day prediction accuracy for the Hubble Space Telescope from July 1996 to 
December 2012. 
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Since late 2012 the JSpOC has been performing the OD for HST using the ASW software. In 
Figure 3, the JSpOC ASW 3-day prediction accuracy for HST is plotted. The 3-day prediction 
accuracy is as high as 8 km in mid-January 2013 and averages 2.35 km over the period from late 
August 2012 to late January 2013. The 27-day rotational period of the Sun is apparent in the 
sinusoidal pattern of the F10.7 index. A key difference between Fermi O/O and HST OD results 
is that GPS point solutions are processed in ODTK for Fermi, while for HST GTDS processed 
ground tracking data and ASW used skin radar tracking. GPS tracking data processed by an 
EKF will likely provide a more accurate definitive orbit solution than ground or skin radar 
tracking since it provides nearly continuous measurement data of higher accuracy. 
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Figure 3: Zoom-in of Hubble 3-day prediction accuracy from JSpOC ASW solutions over 
period from late August 2012 to late January 2013. 


2.2 LEO Drag Modeling 

Beyond the type of measurement data, the other critical difference between the solutions 
generated by the Fermi Operations Team (FOT) and the JSpOC is the drag modeling in the OD 
software. While the FOT uses ODTK to estimate a correction to the atmospheric density, the 
JSpOC ASW software uses what is termed the “High Accuracy Satellite Drag Model” (HASDM) 
[4,5]. The HASDM incorporates Dynamic Calibration of the Atmosphere (DCA) in which 
corrections to the baseline atmospheric density model are estimated via 75-80 calibration 
satellites in varying orbits with the goal of an improved atmospheric density model or OD and 
prediction [6]. Researchers at the University of Kansas have conducted a series of studies 
investigating atmospheric density modeling in ODTK, finding comparable density estimation 
performance to the HASDM with the DCA when precision orbit ephemerides are used [7], This 
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result indicates that even with incorporation of the HASDM, the OD accuracy of the JSpOC 
ASW solutions with low-quality radar tracking is not expected to be better than O/O solutions 
given high-quality GPS measurements. The criticality of drag modeling in ODTK for accurate 
orbit prediction of ICESat, at an altitude that is roughly 60 km higher than Fermi, has also been 
documented by Vallado [8], 

3. Tuning Methodology 

The goal of the orbit determination tuning process for Fermi is to identify the best ODTK filter 
settings, models, and model parameters to improve prediction accuracy. Prediction accuracy is 
closely coupled with definitive accuracy, so the best tuning cases are expected to also be of high 
definitive accuracy. That is, accurate forward predictions generally require a starting state that is 
close to the truth state as the errors will only grow as the spacecraft trajectory is propagated. 
Effectively, the tuning process is an optimization problem where the filter and model parameters 
are the design variables and the prediction accuracy is the objective function to be maximized 
(i.e., minimize error between predictive and definitive states). Three days are selected as the 
prediction span to be “optimized,” as that is roughly when the spacecraft operations team must 
determine whether or not to maneuver to avoid a potential collision with another object. This 
optimization problem is a mixed one, with discrete model types to be selected and continuous 
filter and model parameters to be varied. As most complex optimization problems must be 
solved iteratively, numerical schemes are typically applied. However, the scope of the study 
does not allow for this type of investigation and a manual iterative strategy is employed. A shot- 
gun approach is applied to identify key relationships between the parameter settings and models 
before fine tuning with the objective of improving the 3 -day prediction accuracy over the 
evaluation period. 

A key consideration in evaluating the results is that the tuning was conducted with a two-month 
span of measurement data, from early July 2012 to the end of August 2012. The tuning is, thus, 
inherently biased towards performance over this particular data set and space weather conditions. 
In general, the performance improvement realized with the tuning is not expected to be as 
prominent for all time periods. It is also important to note that the longer the tuning span, the 
more robust the tuning is expected to be (performance -wise) across variable conditions and 
extreme cases in different periods. Wide variability in prediction accuracy is expected over short 
time periods examined (e.g., a few days) given the overall sensitivity of the filter to tuning 
settings. Tuning would likely best be conducted with an even longer period of measurement data 
(e.g., a year instead of two months), but the computational expense of filtering high-frequency 
point solution data is prohibitive for long periods given the iterative nature of the tuning process. 
Furthermore, given the large impact between prediction accuracy during solar maximum and 
solar minimum, adaptive tuning is warranted, where the tuning is updated regularly. 

3.1 Orbit Determination Process Modifications 

A key modification implemented for this study is the use of the ODTK restart capability. In 
ODTK, the filter processing can be stopped (e.g., after the most recent daily measurement data is 
processed) and then started at a later time with the filter initialized to the filter state at the point 
that the restart record was created during the prior run. The current O/O procedure is to start the 
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filter “cold” when it is initialized for the daily OD. That is, the filter covariance information and 
the estimated values for parameters such as ballistic coefficient, atmospheric density correction, 
and solar radiation pressure are not retained after a daily run. The filter starts with only rough 
guesses for the initial covariance matrix and must process measurements to converge back to the 
estimated values. 

For this study, an automated script is implemented to systematically process daily measurements 
in restart mode for successive days. Using the automation script, ODTK filters forward with one 
day of GPS point solution data, creating a restart record at the end of the 24-hour period. The 
filter is then stopped and ODTK is commanded to smooth back over three days of filtered states 
starting from the end of the daily filter run. Finally, starting from the end of the daily filter run 
(or the start of the smoothing process), ODTK’s propagator is used to predict the orbit state 10 
days forward. This automated process generates a daily ephemeris file with three days of 
smoothed definitive states (two days before the filter start and the one day of filter processing) 
and 10 days of predicted states. The process was repeated using a restart record to initialize the 
filter for each day analyzed so that for each day measurements are processed, a separate 
ephemeris file was generated with both standard definitive and predictive states periods. 

In addition to the automated processing with filter restart, other ODTK settings are adjusted to 
follow best practices developed by the Flight Dynamics Facility for Earth orbiting spacecraft in 
ODTK. Namely, the degree and order of the gravity force model was increased from 40x40 to 
70x70. Also, the degree and order of the geopotential used in the variational equations to 
calculate the state error covariance matrix is changed from degree two and order zero to degree 
five and order five. 

3.2 Tuning Strategies 

The tuning process exploits the eleven strategies and variations outlined below. The tuning 
strategies are executed in both combination and isolation to determine the impact of the various 
tactics on OD and prediction accuracy over the two- month evaluation span. Pertinent models 
and parameters that are held constant during the study are listed in Table 1. Models and 
parameters that are varied are detailed in the following Section 3.2.1 - 3.2.1 1. 


Table 1: Fixed ODTK settings for tuning run samples 


ODTK Model/Parameter 

Setting 

Gravity model 

EGM96 

Gravity model degree and order 

70x70 

Degree and order of geopotential variational equations 

5x5 

Spacecraft drag area 

14.18 m 2 

Spacecraft mass 

4357.33 kg 

Drag model 

spherical 
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3.2.1 Measurement Noise 


ODTK enables adjustment of the square root of the variance representing the random Gaussian 
uncertainty of the measurements, or what is termed measurement “white noise sigma.” A white 
noise sigma of 15 meters is used operationally; however, given the high accuracy of GPS point 
solution measurements (ODTK ‘NavSoT measurements) a tighter measurement noise sigma can 
potentially improve definitive accuracy at the cost of filter robustness. To investigate the 
influence of random uncertainty in the measurements on OD accuracy, the white noise sigma was 
varied between 2 and 15 meters for the point solutions. 

3.2.2 Drag Coefficient 

The drag coefficient, Cd, of Fermi is varied between 1.6 and 2.5, where 1.6 is used operationally. 
The drag coefficient is varied to alter the ballistic coefficient value used in the drag model, while 
the drag area is held constant at 14.18 m 2 . In ODTK the ballistic coefficient is defined as 

BC = C D A/m , (1) 

where A is the drag area, and m is the spacecraft mass. Nominally, to calculate the atmospheric 
drag perturbation, the predicted values of the FI 0.7 solar flux and geomagnetic indices are used 
during tuning. However, in two tuning runs, the observed values are used for the atmospheric 
density model so that the drag coefficient value is estimated more accurately. 

3.2.3 Ballistic Coefficient Stochastic Model 

Different stochastic models for the ballistic coefficient are also examined, /?, a parameter 
describing the correction to the nominal ballistic coefficient, is estimated in ODTK, such that 

p = A BC/BC. (2) 

Version 6.1.3 of ODTK offers three built-in stochastic models to model this parameter: Gauss- 
Markov, Vasicek, and random walk. Seago et al. outlines the ODTK implementation of the 
Gauss-Markov and Vasicek stochastic processes [9]. A Gauss-Markov model is used in the 
operational ODTK scenario as it is the only available built-in option in ODTK version 5.0.3 (the 
configured version used currently in the by the Fermi O/O). A Gauss-Markov model is often 
employed for stochastic processes in an OD filter, but the drawback of the model is that the 
estimated value decays back to the nominal value in the absence of measurements, making 
tuning for prediction accuracy more challenging. Later versions of ODTK added a Vasicek 
model, which is composed of both a long- and short-term component. The Vasicek model 
allows for convergence to a nominal mean estimate for the BC (long term component) and then 
corrects on that mean estimate as the filter processing continues (short term component). The 
Vasicek model retains the estimated long-term component through prediction instead of the 
nominal input value. This property makes the model attractive when prediction accuracy is 
important. In the random walk model the estimated value is also held constant by the propagator 
in the absence of measurements, while the uncertainty grows according to a specified diffusion 


8 


coefficient. Johnson compares the three models for prediction accuracy and demonstrates the 
benefit of employing a Vasicek model for the BC using simulated data in ODTK [10]. 

3.2.4 Ballistic Coefficient Uncertainty 

In addition to altering the ballistic coefficient model, the root-variance of the BC correction 
stochastic model is varied. Tightening (i.e., reducing) the a priori uncertainty estimate of BC 
potentially reduces unmodeled atmospheric variations from aliasing into the BC estimate. The 
density variations are instead channeled more appropriately into the atmospheric density 
correction estimate (when both BC and density corrections are estimated simultaneously). 
However, loosening (i.e., increasing) the a priori uncertainty of BC can allow for more 
innovation in the filter. Appropriate variations in combination with the atmospheric density 
correction parameters were investigated. In the Vasicek model, there is both a short- and long- 
term uncertainty, and both of these values are varied when this model is implemented. The BC 
correction sigma is varied between 0.01 and 1.0 in this study. 

3.2.5 Atmospheric Density Model 

Different atmospheric density models are evaluated to determine the best model for Fermi’s 
altitude and the current solar conditions. The baseline model generates a nominal density given 
the predicted solar flux and geomagnetic indices, and a number of model variations each with its 
own advantages and drawback [11]. The Jacchia-Roberts and CIRA 1972 models are evaluated 
in the tuning as they have been shown to be accurate baseline model for low-Earth orbiter OD 

[12] . A correction to the baseline atmospheric density is estimated by the filter in most of the 
tuning runs. However, in several cases, only corrections to the BC are estimated without solving 
for correction to the base density. 

3.2.6 Interpolated versus Step-Function Geomagnetic Indices 

To generate a nominal density, the atmospheric models require a measure of the geomagnetic 
activity. This input is provided as a predicted value of either the a p or K p geomagnetic index 

[13] . K p , the geomagnetic planetary index, is a three-hour, quasi-logarithmic quantity based on 
the mean geomagnetic indices (K indices) from a number of geomagnetic observatories. a p , the 
geomagnetic planetary amplitude, is linear index derived from K p . Predicted values for K p and 
a p are provided in 3-hour steps in the CelesTrak Space Weather file [14], which is updated daily. 
The index used in the atmospheric density model can be an interpolated value based on a cubic 
spline or as a step function input. Both the a p and K p geomagnetic indices are examined as step- 
function and cubic-spline inputs. 

3.2.7 Atmospheric Density and Ballistic Coefficient Correction Half-life Value 

Another critical parameter in the Gauss-Markov and Vasicek stochastic models for the ballistic 
coefficient correction is the exponential half-life value [9]. The half-life defines how quickly the 
state estimate decays back to its nominal value in the absence of measurement data. When 
estimating corrections to the density and BC simultaneously, it is important that the half-life 
values are significantly different so that the two states can be resolved in the filter [15]. These 
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two half-life parameters were critical to prediction accuracy as the atmospheric density is highly 
uncertain and large drag variations from the baseline model are experienced when the solar 
activity is high. Half-life values from 1 minute to 20 days are evaluated. 

3.2.8 Atmospheric Density Correction Sigma Scale 

An experimental setting available to the user in ODTK is a multiplicative factor applied to the 
uncertainty of the atmospheric density correction. This parameter called the “density correction 
sigma scale” is applied at all times during fdtering, and can be increased to produce a looser 
filter with the potential for greater innovations. During periods of high solar activity, a 
sufficiently loose filter can better estimate large density fluctuations that cannot easily be 
predicted and, thus, are not present in the model. The risk of increasing the sigma scale is that 
the definitive accuracy is adversely affected. The density sigma scale is varied between 0. 1 and 
4.0 in this investigation. 


3.2.9 Atmospheric Density Ratio Increase Threshold 

As described in the ODTK Help File [16], the atmospheric density ratio increase threshold: 

specifies a threshold for the stepwise change of a ratio of the density evaluated at perigee 
at the current solar and geomagnetic conditions to an evaluation of the density at mean 
solar and geomagnetic conditions. When this threshold is exceeded in the filter, the 
process noise for atmospheric density changes from the baseline model to the active sun 
model to open up the covariance under conditions of high solar activity. 

The density ratio increase threshold is tested at values of 0. 1 and 1 .0. In the best tuned cases, the 
mean prediction accuracy is not strongly affected by this parameter. However, in certain tunings 
this threshold has a notable influence on OD and prediction performance. Specifically, when it is 
coupled with the type of geomagnetic index employed. 

3.2.10 Data Thinning 

Operationally, the Fermi Operations Team receives GPS point solution measurements at a one- 
second frequency. To reduce the time required to perform OD, the data used within the 
operational filter is down sampled to one point solution per minute. Thinning the data, however, 
not only affects the computational expense of the filter, but also the fundamental state of the 
filter at each step due to the stochastic models representing key parameters. By thinning the 
data, the estimated value of the ballistic coefficient correction and atmospheric density correction 
can be retained for a longer time through the prediction with appropriate stochastic model 
parameter settings. This approach enables long exponential half-life values to be used while still 
allowing for filter flexibility. ODTK allows for variable thinning and 1 sec, 1 min, 5 min, 8 min, 
10 min, 12 min, and 20 min thinning rates are examined. 
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3.2.11 Integrators 


The final tuning strategy is the choice of variable-step integrators. A Runge-Kutta-Fehlberg 
integration method of 7 th order with 8 th order error control is evaluated against a Runge-Kutta- 
Vemer method of 8 th order with 9 th order error control for the integration step size. Additionally, 
the minimum step size of the integrators is tested at 0.5 and 1 seconds, and the relative error 
tolerance is varied as well. 

3.3 Tuning Process and Insights 

The numerous tuning options outlined in Section 3.2 are investigated in detail, where the results 
of each tuning run inform the subsequent runs. The definitive position and velocity uncertainties 
are relatively unaffected regardless of the tuning. The position uncertainty does increase when 
the data is thinned, but it is not suspected that there is a significant change in the definitive 
accuracy given the density and continuity of the point solution measurements. While there is no 
truth data to compare against, this supposition is supported by the fact that thinning the data 
often improved prediction accuracy. The problem is then less of an orbit determination problem 
and more of an orbit prediction problem. This realization supports a tuning approach with the 
goal of improving average prediction accuracy over the evaluation period. A secondary goal of 
minimizing the maximum 3 -day error is also considered, but this objective is not always in 
competition with the primary objective of minimizing the mean 3 -day position error over the two 
month evaluation period. That is, the mean 3 -day error and maximum 3 -day error frequently 
increased or decreased in step relative to the other runs. 

Over 130 ODTK filter/prediction runs are conducted in this study using a two month evaluation 
span from early July 2012 to the end of August 2012. Widely varying results in 3 -day position 
prediction accuracy are observed as illustrated in Figure 4. The earlier tuning runs are in dark 
blue and the later runs are in light blue. The best tuning in terms of average prediction error is 
highlighted in red. Note that plotted in green on the secondary axis is the exospheric 
temperature. This index is calculated from the FI 0.7 flux and the geomagnetic indices in the 
daily space-weather file that ODTK ingests. The exospheric temperature is input into the 
atmospheric model (Jacchia-Roberts or CIRA 1972) to generate the density used in the drag 
calculation. 
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Figure 4: 3-day prediction accuracy for all tuning runs over the two month evaluation 
period. Also plotted on the secondary axis is the exospheric temperature in Kelvin. 


The correlation between exospheric temperature and prediction accuracy is apparent in Figure 4. 
The exospheric temperature provided a gauge of atmospheric density variations and is a 
reasonable predictor of 3 -day prediction accuracy. That is, when the exospheric temperature is 
high or there were dramatic fluctuations, there is generally a more significant mismodeling of the 
density through the predictions resulting in higher prediction errors. Additionally, the later runs 
in light blue (and the best-tuned run in red) more closely track the exospheric temperature, 
especially in the last week of August when the difference between the worst and best tunings is 
largest. Note that some of the earlier runs in dark blue exhibit better predictive accuracy during 
certain spans of the two month evaluation period (e.g., from 7/26/12 - 7/29/12 and from 8/9/12 - 
8/19/12), but produce relatively poor mean accuracy over the entire period due to large 3 -day 
prediction errors during other times (e.g., from 7/30/12 - 8/8/12 and from 8/23/12 - 8/30/12). 


The effect of the tuning strategies is most prominent when parameters and models are varied in 
combination. Plotted in Figure 5 is the mean 3 -day prediction accuracy over the evaluation 
period as a function of tuning run, where the run index was increased with each run 
chronologically. Early runs in the tuning analysis focus on exploration of a diversity of settings 
as is apparent in the early fluctuations of mean accuracy, while the later runs focus on refinement 
and exploitation. Runs 78 and 79, the two cases with the best mean accuracy, employ observed 
values for the FI 0.7 solar flux and geomagnetic indices to help refine the nominal drag 
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coefficient (predicted values are used for all other runs). From this result, better accuracy is 
achieved with an improved atmospheric model, one that more closely represents the environment 
in which the vehicle is flown. However, even with an improved representation of the 
atmosphere, the 3-day prediction errors are on the order of 300 m, or about 25% lower than the 
best tuning with predicted values. 
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Figure 5: 3-day prediction accuracy as a function of tuning run. 


4. Tuning Results 

This section highlights several tuning runs, detailing some of the findings and demonstrating 
how the prediction accuracy is affected by the parameter settings. Five different representative 
tunings are discussed (runs 91, 32, 50, 68, and 117), and the associated critical ODTK settings 
are outlined in Table 2 (N/A in the table indicates that the parameter is not applicable). Figure 6 
illustrates the daily 3 -day prediction plotted over the two month evaluation period along with the 
exospheric temperature indicated by the red line for the five cases. 
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Table 2: < 

3DTK settings for tuning run samples 

Parameter 

ODTK Setting 
Run 91 (O/O) 

ODTK Setting 
Run 32 

ODTK Setting 
Run 50 

ODTK Setting 
Run 68 

ODTK Setting 
Run 117 

Data frequency 

1 min 

10 min 

10 min 

10 min 

10 min 

BC model 

Gauss- 

Markov 

Vasicek 

Vasicek 

Random walk 

Gauss 

Markov 

BC half-life 

7 days 

2 days 

7 days 

N/A 

10 min 

BC short-term a 

0.3 

0.5 

1.0 

0.05 

0.1 

BC long-term a 

N/A 

0.5 

1.0 

N/A 

N/A 
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Figure 6: 3-day prediction accuracy for runs 117, 68, 50, 32 & 91 as a function of epoch. 
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4.1 0/0 ODTK Settings (run 91) 


To provide a benchmark for tuning comparison, the settings from the owner/operator operational 
ODTK scenario are used in run 91. The most critical settings from run 91 are listed in Table 2. 
A fundamental difference between the actual operational results and the results of run 91 is that 
the Fermi flight operations team (FOT) does not the follow the same OD procedure used in this 
analysis as discussed in Section 3.1. Additionally, the FOT uses ODTK 5.0.3, which has fewer 
tuning options than ODTK 6.1.3 used in this analysis. The key limitation between the versions 
is that only a Gauss-Markov ballistic coefficient is available in ODTK 5.0.3. Operationally, the 
GPS measurements are available at one measurement per second, but the data is thinned to one 
measurement per minute. The mean 3-day prediction accuracy from run 91 is 624 m. The worst 
daily prediction is from 7/15/12, which correlates to the highest exospheric temperature during 
the evaluation period. Additionally, relatively poor accuracy is achieved at the end of August 
when the exospheric temperature is relatively low as illustrated in Figure 6. Given the expected 
correlation of prediction accuracy with the exospheric temperature, accuracy may be improved 
over this span. 

The ballistic coefficient correction, the estimated change divided by nominal input value, and 
atmospheric density correction for run 91 appear in Figure 7 and Figure 8, respectively. Recall 
from Table 2 that the half-life for the BC correction is relatively long compared to the 
atmospheric density half-life (7 days compared to 180 minutes). The large disparity between the 
two values is critical to separating the two channels in the filter. However, the BC varies by up 
to 100% at the end of the two-month period and the secular, long-term variation is much more 
prominent than any orbit-to-orbit changes as shown in Figure 7. The atmospheric density 
correction is much more muted with short terms variations as indicated in Figure 8. This 
behavior is not expected, as in reality the BC likely only varies orbit-to-orbit with little secular 
drift, while the atmospheric density correction is expected to fluctuate more widely and with 
longer term deviations due to the changing space weather that strongly affects the density. A 
tuning that enables the filter to estimate corrections that better reflect the actual physics is 

desired as it will likely be more robust. 

Ballistic Coefficient Correction 



Ballistic Coefficient Correction dB/B dB/B +2 Sigma 

dB/B -2 Sigma 

Figure 7: Ballistic coefficient correction as a function of epoch for run 91. 
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Atmosphere Density Correction 



Atmosphere Density Correction dp/p dp/p +2 Sigma 

dp/p -2 Sigma 

Figure 8: Atmospheric density correction as a function of epoch for run 91. 

4.2 Thinning Data (run 32) 

One of the key findings from the tuning analysis is that thinning the data can be very effective. 
In run 32, the data was thinned down to one measurement every 10 minutes. This longer period 
between measurements permits lengthening the half-life values for ballistic coefficient or density 
correction while allowing for the filter to incorporate innovation in the estimation. 

In addition to thinning the measurement data, different stochastic process models for the BC 
correction are evaluated. A Yasicek, not Gauss-Markov, model is used in run 32. The Vasicek 
stochastic process enables estimation of both a short-term and long-term bias from an a priori 
BC. It is appealing for the BC correction model because the long-term time -varying bias 
estimate does not decay in the absence of data, while the short-term bias estimate decays at a rate 
dictated by the associated half-life value. This process allows for the long-term offset from the a 
priori to be carried through the prediction, while still allowing for the estimation of orbit-to-orbit 
variations from the estimated long-term value via the short-term component of the model. In run 
32 the BC correction uncertainty is also increased, while a shorter density correction half-life 
value of 360 min is tested. Additionally, the Cd is increased substantially from 1.6 in run 91 to 
2.4. The mean 3-day error decreases from 624 m to 498 m primarily because of the introduction 
of data thinning and increased BC uncertainty. However, as in run 91, the BC estimation 
absorbs the majority of the drag variation, which is likely not physical. The estimated 
corrections to BC and atmospheric density for run 32 are shown in Figure 9 and Figure 10, 
respectively. 
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Ballistic Coefficient Correction 



Ballistic Coefficient Correction dB/B dB/B +2 Sigma 

dB/B -2 Sigma 

Figure 9: Ballistic coefficient correction as a function of epoch for run 32. 
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Figure 10: Atmospheric density correction as a function of epoch for run 32. 


4.3 Stochastic Model Parameter Testing (run 50) 

In ran 50, the BC correction half-life is extended to seven days, as compared to two days in ran 
32. This adjustment enables the estimated BC correction to be carried through a longer 
prediction period before the BC decays to long-term bias estimate. The atmospheric density 
half-life is also lengthened to one day for the same purpose. Additionally, the BC correction 
sigmas in the Vasicek stochastic model are increased to open the covariance of the filter. 
Marked improvement in performance in the last part of August is apparent in comparing ran 50 
(orange) to runs 91 (light blue) and 32 (green), where the mean 3-day prediction error is reduced 


17 


to 436 m. As illustrated in Figure 11 and Figure 12, the BC correction estimate effectively 

absorbs the majority of the drag correction, which is not desired. 

Ballistic Coefficient Correction 



Ballistic Coefficient Correction dB/B dB/B +2 Sigma 

dB/B -2 Sigma 

Figure 11: Ballistic coefficient correction as a function of epoch for run 50. 
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Figure 12: Atmospheric density correction as a function of epoch for run 50. 


4.4 Stochastic Model Variation and Half Life Value Reversal (run 68) 

A random walk stochastic model for BC is employed in run 68. This type of stochastic process 
allows for the estimated value of BC to be strictly maintained without decaying through the 
prediction. This model is beneficial in terms of prediction accuracy as shown by the dark blue 
line in Figure 6. Another significant change in this run is extending the atmospheric density 
correction half-life to ten days. In using a long density correction half-life with a random walk 
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BC, it is clear that the density correction estimate varies widely instead of the BC as depicted in 
Figure 13 and Figure 14. Corrections to the BC are estimated, but the density correction absorbs 
the majority of the drag correction as desired. To encourage innovation in the density correction 
channel, and minimize the potential for filter smugness, the density sigma scale is increased from 
the nominal value of 1.0 to 3.0. The inflation of atmospheric density uncertainty helps to guide 
any density variations from aliasing into the BC correction estimation. This tuning shows further 
improvement in the prediction accuracy at the end of August 2012, and more closely follows the 
exospheric temperature variations. The mean 3-day prediction error is further reduced to 412 m 
with this tuning strategy. 


Ballistic Coefficient Correction 



Ballistic Coefficient Correction dB/B dB/B +2 Sigma 

dB/B -2 Sigma 

Figure 13: Ballistic coefficient correction as a function of epoch for run 68. 


Atmosphere Density Correction 



Atmosphere Density Correction dp/p dp/p +2 Sigma 

dp/p -2 Sigma 

Figure 14: Atmospheric density correction as a function of epoch for run 68. 
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4.5 Refinement (run 117) 


Given that the ODTK version used by the owner/operator only offers a Gauss-Markov stochastic 
model for BC, further tuning focused on identifying appropriate parameter settings for a Gauss- 
Markov BC stochastic process. A Gauss-Markov stochastic model is employed in run 117 with 
the goal of tuning parameters to achieve an average 3-day prediction error in line with run 68. A 
Gauss-Markov BC model, however, does not retain its estimated value through prediction and 
instead decays towards the a priori value according to the specified half-life without any 
measurements. The parameters for run 117, the best tuned case in study, are listed in Table 2 
and the 3-day prediction accuracy is compared against the other runs described illustrated in 
Figure 6. This tuning is the result of extensive refinement of all the tuning options identified in 
Section 3.2. A short BC correction half-life of 10 minutes and long atmospheric density 
correction half-life of 8 days guide the drag correction variations into the density correction 
channel as illustrated in Figure 15 and Figure 16. The short BC correction half-life does not 
allow for significant correction of the BC as the data frequency is also thinned to 10 minutes. 
The BC correction quickly decays towards the input value by the time of the next measurement, 
hindering the filter from estimating large corrections to the BC. The scale of Figure 15 does not 
indicate that there was any correction to the BC; however, corrections on the order of .1% do 
exist. Only minor corrections to the BC necessitate that the initial BC value is appropriate for 
predicting forward, i.e., close to actual spacecraft BC. The initial BC is in this run is based on 
the BC estimated in run 78 when the observed solar flux and geomagnetic indices are used. 

The mean 3-day position error of 41 1 m characterized in Figure 6 is slightly better than the 412 
m performance of run 68. The maximum error in run 117 on July 15 is also slightly lower than 
that of run 68. Interestingly, there is a five day period (8/12/12 - 8/16/12) where the 3 -day 
prediction accuracy of runs 91 and 32 are better than run 117. This behavior demonstrates that 
the best tunings are only better on average, and also indicates that tuning over long periods 
allows for an OD system that will provide robust performance over a broad set of conditions. 
Adaptive tuning, however, is expected to provide the best performance, where regular tuning 
updates are better adjusted to the current operating and environmental conditions. 
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Figure 15: Ballistic coefficient correction as a function of epoch for run 117. 
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Atmosphere Density Correction 



Atmosphere Density Correction dp/p dp/p +2 Sigma 

dp/p -2 Sigma 

Figure 16: Atmospheric density correction as a function of epoch for run 117. 

5. Three Month Results 

To further compare the best tuning against the 0/0 ODTK settings, ODTK runs using the 0/0 
settings (run 91) and the best tuned case (run 117) are examined over the approximately three 
month period for which Fermi measurement data is available. The 3-day prediction accuracy 
comparison from 7/6/12 to 9/30/12 appears in Figure 17. 


3000 



Solution Epoch 


Figure 17: Comparison of 3-day accuracy between run 91 (0/0 settings) and run 117 (best 
tuned case) from July 6 to September 28 2012. 
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In the first part of July 2012, the performance of the 0/0 settings of run 91 does not dramatically 
differ from run 117, the best tuned cased. However, there are periods where the tuned ODTK 
process performs much better than the 0/0 settings, notably from late July into early August, 
during the end of August, and again at the end of September. Significant performance 
improvements are realized, where a 45% improvement in the mean 3-day accuracy is achieved as 
well as a reduction in the maximum error by 26%. The tuning is especially effective in 
improving the prediction accuracy on the worst-case days from run 91 tuning, with the exception 
of July 15. A statistical comparison of the 3-day accuracy between run 91 and run 117 is 
su mm arized in Table 3. The 26% improvement in standard deviation of the 3-day prediction 
accuracy over the three months indicates that not only are the predictions more accurate on 
average, but they are also more consistent. 

Table 3: Statistical comparison of 3-day accuracy between run 91 and run 117 from July 6 


to September 28 

!, 2012 


0/0 ODTK Settings 

Refined ODTK Settings 
Run 117 

Improvement 

Mean (m) 

755 

425 

45% 

Maximum (m) 

2923 

2164 

26% 

Standard deviation (m) 

609 

390 

36% 


Significant improvement in longer-term predictive accuracy is also realized, where the run 117 
tuning strategy reduces the 6-day mean error by 36% over 0/0 settings as illustrated in Figure 18 
and Table 4. Similar to the the 3-day prediction accuracy, the most substantial improvement is 
during the end of August and the end of September. The run 117 tuning also reduces the 
maximum 6-day prediction error on September 1, 2012 from 14.8 km to 9.6 km, a meaningful 
improvement. Gains are also attained in prediction consistency, where the standard deviation of 
the 6-day prediction accuracy over the three months is reduced by 42%. 
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Figure 18: Comparison of 6-day accuracy between run 91 (0/0 settings) and run 117 (best 
tuned case) from July 6 to September 25 2012. 
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Table 4: Statistical comparison of 6-day accuracy between run 91 and run 117 from July 6 


to September 28 2012 



0/0 ODTK Settings 

Refined ODTK Settings 
Run 117 

Improvement 

Mean (m) 

3412 

2161 

37% 

Maximum (m) 

14809 

9571 

35% 

Standard deviation (m) 

3007 

1751 

42% 


The 1-day prediction accuracy between runs 91 and 117 is also examined over the 3 -month 
period. A 50% error reduction in mean 1-day error is achieved, larger than the improvement in 
the 3- and 6-day accuracy. It is apparent that the predictions are more consistent in run 117 as 
illustrated in Figure 19. This consistency is reflected in the reduced standard deviation of the 1- 
day prediction accuracy listed on Table 5. 



Figure 19: Comparison of 1-day accuracy between run 91 (0/0 settings) and run 117 (best 
tuned case) from July 6 to September 25 2012. 


Table 5: Statistical comparison of 1-day accuracy between run 91 and run 117 from July 6 
to Se ptember 28 2012 



0/0 ODTK Settings 

Refined ODTK Settings 
Run 117 

Improvement 

Mean (m) 

87 

43 

50% 

Maximum (m) 

230 

145 

37% 

Standard deviation (m) 

61 

32 

46% 
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6. Formal Error Envelope 


The 3-day position prediction error is also evaluated on a component-by-component basis. 
Representative radial, in-track, and cross-track 3 -day position errors over the two month 
evaluation period are plotted in Figure 20, Figure 21, and Figure 22 respectively. In addition to 
the prediction error for each position component, the 3 -sigma formal error envelopes derived 
from the propagated covariance matrix after a 3 -day prediction are plotted in black. As 
expected, the in-track component is the largest contributor to the position error. The in-track 
error also corresponds to the exospheric temperature metric as illustrated in Figure 21. The in- 
track error remains conservatively within the formal 3 a envelope even when there are large 
jumps in exospheric temperature. The radial and cross-track component errors are significantly 
smaller, and are also bounded by the 3 -sigma formal error envelope. 
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Figure 20: Representative 3-day prediction radial position error as function of epoch 
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Figure 21: Representative 3-day prediction in-track position error as function of epoch 
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Figure 22: Representative 3-day prediction cross-track position error as function of epoch 


Figures 20 through 22 highlight possible biased components. Neither the radial nor the in-track 
3 -day errors reveal a bias on either side of zero, but a roughly seven-day bias is evident in the 
cross-track error shown in Figure 22. This bias is likely due to stagnant values in the daily 
predictive Earth Orientations Parameter (EOP) file. The predicted values in the daily EOP file 
are generated from two sources with different predictions and irregular updates. In the EOP files 
evaluated over this 2 month span, there were several days of stagnant predict values. The same 
prediction values over multiple days are the likely source of the cross-track error bias. However, 
a critical point is that these biases have very little effect on the predictive accuracy. 


7. Conclusion 


A tuning analysis to improve the orbit determination and prediction accuracy of Fermi reveals 
that the differences between Fermi owner/operator and JSpOC orbit solutions is primarily a 
prediction issue rather than an orbit determination issue. Given the relatively large effect of 
atmospheric drag on the spacecraft, it is difficult to predict the orbit multiple days in to the future 
with high accuracy because of the uncertainty and variability of the atmospheric density in the 
drag model. The definitive position uncertainty does not vary widely and higher definitive state 
uncertainties do not correlate with poorer prediction performance. Additionally, the prediction 
accuracy of Fermi is consistent with that of satellites with a similar altitude. For the dates 
examined, the Fermi 3-day prediction accuracy using a variety of tools is better on average than 
those of Hubble Space Telescope, which resides at a higher altitude. While the prediction 
performance of the O/O solutions examined are not overtly poor relative to other satellites at 
similar altitudes, substantial improvements can be attained with different ODTK tunings. An 
extensive tuning analysis produced an average 3 -day prediction error that is 45% lower than the 
tuning in the current O/O ODTK scenario over a 3-month period in 2012. 
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Further improvements in prediction accuracy may feasible. Potential future work to reduce the 
prediction error beyond the tuning analysis in this study includes incorporating attitude- 
dependent drag modeling in the filter (e.g., box-wing model). Additionally, OD definitive 
accuracy might be improved with application of a GPS-elevation constraint to reduce ionosphere 
modeling error. Further improvement in definitive accuracy may possible by including one or 
more additional satellites, and their associated measurement data, in the filter to better resolve 
atmospheric density variation. Any improvement in definitive accuracy is likely to translate into 
improved prediction accuracy. 
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